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Abstract 

Transient grating spectroscopy has emerged as a useful technique to study thermal phonon 
transport because of its ability to perform thermal measurements over length scales comparable to 
phonon mean free paths (MFPs). While several prior works have performed theoretical studies of 
quasiballistic heat conduction in transient grating, the analysis methods are either restricted to one 
spatial dimension or require phenomenological fitting parameters. Here, we analyze quasiballistic 
transport in a two-dimensional transient grating experiment, in which heat conduction can occur 
both in- and cross-plane, using an analytic Green’s function of the Boltzmann equation we recently 
reported that is free of fitting parameters. We demonstrate a method by which phonon MFPs can 
be extracted from these measurements, thereby extending the MFP spectroscopy technique using 
transient grating to opaque bulk materials. 
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I. INTRODUCTION 


Transient grating spectroscopy (TG) is a noncontact, optical pump-probe experiment that 
is capable of measuring the in-plane thermal conductivity of thin films and bulk materials. 1,2 
In this experiment, the sample is impulsively heated with a spatially sinusoidal temperature 
profile created by splitting a pump laser pulse using diffractive optics and recombining the 
beams on the sample. A probe laser subsequently observes the in-plane transient thermal 
decay by measuring the intensity of light diffracted from the sample surface. In the heat 
diffusion regime, the shape of the transient decay curve is controlled by the in-plane thermal 
diffusivity, 2 and thus the thermal diffusivity can be measured from the exponential time 
constant of the thermal decay. 

Recently, considerable interest has focused on quasiballistic heat conduction in TG, in 
which phonon mean free paths (MFPs) are comparable to the grating period. 3-8 Prior works 
have demonstrated that observations of thermal transport in this regime yield information 
about the MFP spectrum of thermal phonons in a technique known as MFP spectroscopy. 5 ’' 
A number of models have been proposed to explain these results, including a two-fluid 
model, 5 a two-parameter nondiffusive model, 3 and related models for time-domain thermore¬ 
flectance experiments. 9-12 Recently, we reported an exact analytic solution of the Boltzmann 
transport equation (BTE) for ID transport in TG. 8 

However, often the thermal transport induced by the pump pulse is not purely ID due to 
the finite optical penetration length of light in the sample. 2 For example, applying TG to a 
GaAs wafer generates the externally defined in-plane grating but also an exponential decay 
in the cross-plane direction, resulting in both cross-plane and in-plane heat conduction and 
a thermal decay that is no longer a pure exponential. 13,14 While in the diffusive case these 
measurements can be readily interpreted with Fourier’s law, 13 in the quasiballistic regime 
interpreting the thermal decays is difficult because a unique thermal time constant can no 
longer be clearly identified. 

In this work, we provide a theoretical framework for treating multidimensional heat con¬ 
duction in TG based on an analytic Green’s function for the BTE we recently reported. 15 
We End that the transport in 2D TG is considerably more complicated than that in ID 
because the heating profile contains a spectrum of spatial frequencies rather than a single 
spatial frequency as in the ID case. Nevertheless, we demonstrate a method to recover the 
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MFP spectrum from these measurements without any fitting parameters or knowledge of 
the MFPs. This work extends the capability of MFP spectroscopy using TG from thin films 
to bulk opaque solids. 


II. THEORY 


A. Boltzmann Equation 


We begin by briefly reviewing the derivation of the multidimensional Green’s fuction to 
the BTE as reported by Hua and Minnich. 15 The BTE under the relaxation time approxi¬ 
mation is given by: 


dg k 

dt 


+ v ■ V r g k 


g*-go(T(r,t)) | 

Tk 


( 1 ) 


where g k = hcu(f k (r,t,k.) — fo(T 0 )) is the desired deviational distribution function, Q k (r,t) 
is the spectral volumetric heat generation, v is the phonon group velocity, r k is the phonon 
relaxation time, and k is the phonon wavevector in phase space. go(T(r,t)) is the local 
equilibrium deviational distribution function and is given by: 


go{T{r,t)) = hiu(f BE (T(r,t)) - f 0 (T 0 )) « C k AT(r,t). (2) 


assuming a small temperature rise AT(r, t) = T(r, t)—T 0 relative to a reference temperature 
To. Here, h is the reduced Planck constant, cu(k) is the phonon frequency given by the 
dispersion relation, f B E(T(r,t )) is the local Bose-Einstein distribution, / 0 (T 0 ) is the Bose- 
Einstein distribution at T 0 , and C k = k B (x/ sinh(x)) 2 is the mode specific heat, where 
X = hu(k) /2k B T 0 . To close the problem, energy conservation is used to relate g k to AT(r, t) 
as 
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where the sum is performed over all phonon modes in the Brillouin zone. 

An analytic solution to the above equations is obtained by performing a Fourier trans¬ 
form over space and time in Eq. 1, substituting the result into Eq. 3, and solving for the 
temperature response AT(r), £ y , £ z ), where r/ and £ x , £ y , £ z are the Fourier variables in the 
time and spatial domains, respectively. This procedure assumes that the thermal transport 
occurs in an infinite domain so that a Fourier transform can be performed. The result for 
an isotropic crystal is: 
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where the magnitude of spatial frequency £ = \/£|~+£|~+£j. In the weakly quasiballistic 
regime where tjt^ <C 1, the result can be simplified by expanding all terms in powers of r)T k 
and keeping only linear terms. The result is: 
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where Q(rj, £ x , £ x , £ z ) is the heating profile summed over all phonon wavevectors and /c(£) 
is the apparent thermal conductivity that depends only on £, reflecting the isotropy of the 
crystal. The apparent thermal conductivity is given by: 
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where /(A) is the MFP distribution, A is the MFP, S(Kn ) is the suppression function, 
and the Knudsen number Kn = £A. 

We now specialize this result for the 2D TG experiment. Along the in-plane, or x direc¬ 
tion, the heating profile is a sinusoid with spatial frequency £ x . Along the cross-plane, or y 
direction, the heating profile is a decaying exponential with decay length /3\ corresponding 
to the optical depth for the pump beam in the solid. The probe beam samples this temper¬ 
ature profile with an optical depth /3 2 corresponding to the probe wavelength. The heating 
along the z direction is uniform and hence £ z = 0. Recalling the expression for the Fourier 
transform of a decaying exponential, the temperature profile in 2D TG can thus be written 
as: 


A T( V ,Q = A 

where A is an arbitrary constant, 
domain as: 
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This equation can be easily converted to the time 
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where a(£) = k(£)/C is the apparent thermal diffusivity and A' = A/C. 

Equation 9 shows that the transient temperature profile for a 2D grating is a multiex¬ 
ponential decay with a spectrum of time constants determined by the apparent thermal 
conductivities. The relative contribution of each exponential decay to the total signal is 
weighted by the Fourier transform of the cross-plane heating profile. This situation is con¬ 
siderably more complicated than in the ID case, for which the sole spatial frequency in 
the problem is externally defined by the grating period. In the ID case, a single apparent 
thermal conductivity can be identified from the exponential time constant of the observed 
thermal decay for a specific grating period, and a set of thermal conductivities versus grating 
period can then be used to directly recover the MFP spectrum. 7 

In this present 2D case, however, the thermal decay is no longer a single exponential 
but a multiexponential. Such multiexponential decays occur often in science, for example in 
time-resolved measurements of fluorescent decays, 16 and are challenging to interpret because 
recovering all the exponential time constants from a multiexponential decay is an ill-posed 
problem for which a unique solution does not exist. Thus in the 2D case, an inverse problem 
must be solved to recover the apparent thermal conductivities from the measured temper¬ 
ature decays before any further information about the thermal phonons can be extracted. 
Further, the inverse problem is nonlinear because the function to be recovered is inside an 
exponential. 


B. Optimization Problem 

This inverse problem can be solved using constrained nonlinear optimization methods, 
thereby yielding the apparent thermal conductivities as a function of spatial frequency £. 
These values can then be directly used as input to the original inverse problem of recovering 
the MFP spectrum described in Ref.' that is readily solved using convex optimization. 

The present inverse problem is to recover the function a(£) from a set of transient ther¬ 
mal decays, A T(t,£ x ), where the spatial frequency along the in-plane direction can be 
externally controlled by varying the grating period. The optical decay lengths f3± and are 
assumed to be fixed by the material and optical wavelengths used in the experiment, thereby 
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specifying the spatial frequencies along the cross-plane direction, £ y . This prototypical prob¬ 
lem has been addressed using a number of methods in the past, including Prony’s method, 17 
Fourier methods, 18 and linearization to enable linear optimization algorithms. 16,19,20 In this 
work, we use nonlinear optimization methods to solve the problem. 

Let us discretize this problem by assuming with have N transient data sets A Ti(t,£ X j) 
corresponding to different grating periods used in the experiment. Our goal is to adjust M 
apparent thermal diffusivities ay that depend on spatial frequency to match Eq. 9 to the 
experimental data set. This problem can be classified as a nonlinear constrained optimization 
by defining a cost function that measures the difference between the simulated and actual 
thermal decays and attempting to minimize this cost function by varying ay. As with many 
inverse problems, attempts to solve this problem without any constraints will likely not 
succeed because a unique solution does not exist. The desired solution can only be recovered 
from the measurements if known physical facts about the solution are incorporated into the 
optimization. 

Fortunately, a number of facts are known about the solution. Based on the behavior of 
the derivatives of the suppression function, one can derive that the first derivative of a(£) 
must be negative, the second derivative must be positive, and the third derivative must 
be negative. These constraints can be implemented by approximating the derivative with 
finite differences and incorporating the resulting equations into the optimization as linear 
inequality constraints. 

Second, the solution cannot have jumps between adjacent ay that are too large. Consid¬ 
ering that the spatial frequencies are uniformly spaced, the maximum difference between 
oq+i and ay is: 


/ °o / dS \ 

a(A) f -~jj^ ) AdAA£ < a bu i k A med A£ (10) 

where a (A) = k(A)/C is the differential thermal diffusivity distribution versus MFP and 
A med is an estimate of the median MFP of the material. This constraint can also be incor¬ 
porated into the optimization as a linear inequality constraint. 

Third, the apparent thermal diffusivities do not vary substantially over the range of 
experimentally attainable grating periods, and thus upper and lower bounds on the allowed 
values of ay can be imposed. For instance, for the analysis performed for Si in the next 
section, the apparent thermal diffusivities must be between 1 and 100 mm 2 /s. Similar 
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FIG. 1: (a) Synthesized temperature decay curves (solid lines) and decay curves corresponding 
to the thermal diffusivity function reconstructed by nonlinear optimization (dashed lines). The 
reconstructed thermal decays are in good agreement with the original synthesized curves, (b) 
Reconstructed thermal diffusivity a(£) (symbols) with equality constraints compared to the actual 
function (solid line). The reconstructed function is in good agreement with the actual function, 
(c) Reconstructed MFP spectrum using the thermal diffusivity result from (b). The reconstructed 
MFP spectrum and the actual distribution agree well with each other. 


bounds can be imposed for other solids for which basic knowledge of bulk thermal properties 
is known. 

With these constraints, the optimization problem can be solved with standard optimiza¬ 
tion routines in MATLAB. We use the fmincon function in MATLAB for this work. The 
cost function to be minimized is defined as the magnitude of the difference between the 
calculated and measured temperature decays for N different grating periods £ x . The con¬ 
straints described above are put into matrix form and included in the optimization as linear 
inequality constraints. The outputs of the optimization are M discrete values of the function 
a(£j) that minimize the cost function while satisfying the constraints. 


III. RESULTS 

To demonstrate the reconstruction, we first synthesize temperature decay curves AT* for 
different grating periods A * = 2n/£ x ,i between 500 nm to 20 microns. The decay curves 
are obtained by calculating the apparent thermal diffusivities versus spatial frequency using 
Eq. 6 and the phonon dispersion and lifetimes for Si at room temperature as calculated 
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by Jesus Carrete and Natalio Mingo using ShengBTE 21,22 and Phonopy 23 from interatomic 
force constants obtained with VASP. 24 ~ 2 ' These results are then inserted into Eq. 9, yielding 
transient decay curves for a particular grating period. 

The resulting decay curves, with random noise added to simulate experimental uncer¬ 
tainty, are shown in Fig. la. We then use the optimization procedure described in the 
previous section to reconstruct the apparent thermal diffusivities. For this problem, we dis¬ 
cretized the spatial frequencies into M — 30 equally spaced points, and attempted to fit 
IV = 25 temperature decay curves. The extinction coefficients in the cross-plane direction 
are taken to be /3i — fa — 1/500 nm _1 . 

In the first set of results, we took the minimum grating period to be 500 nm, a smaller 
value than is typically achieved experimentally but still possible with commercially available 
optics. The availability of data at a grating period that is comparable to the optical decay 
length is very helpful because the thermal decay corresponds predominantly to the single 
spatial frequency defined by the imposed grating rather than a spectrum of spatial frequen¬ 
cies as in the general case. As shown in Fig. lb, the reconstructed thermal diffusivity is in 
excellent agreement with the actual thermal diffusivity despite the presence of noise in the 
original temperature decay curves in Fig. la. 

Using this reconstructed thermal diffusivity function, we then close the problem by per¬ 
forming the convex optimization procedure described in Ref. 7 to recover the desired MFP 
spectrum. This result is given in Fig. lc, and shows that the reconstructed distribution 
is in good agreement with the actual MFP spectrum. Our approach of solving two inverse 
problems, starting from the original transient decay curves, is thus successfully able to recon¬ 
struct the desired MFP spectrum without any knowledge of the MFPs or thermal properties 
of the material. 

The reconstruction without data at grating periods comparable to the optical decay 
length is more challenging but still feasible. The reconstructed thermal diffusivity versus 
spatial frequency obtained when only using grating periods from 1-20 microns is shown in 
Fig. 2a. The reconstruction and correct answer are in generally good agreement, though 
a discrepancy exists at high spatial frequencies. Performing the second inverse problem 
with this input results in a poor reconstruction of the MFP spectrum for values of the 
MFP below 1 micron (not shown). To circumvent this limitation, we recognize that the 
reconstructed values for spatial frequencies larger than that corresponding to the minimum 




FIG. 2: (a) Reconstructed thermal diffusivity with data from grating periods between 1-20 microns 
(symbols) compared to the actual function (solid line). The reconstruction is good but some 
discrepancy exists at higher spatial frequencies, (b) Reconstructed MFP spectrum (green circles) 
using the only the thermal diffusivities in (a) with spatial frequency less than 27r /im _1 (red 
squares) compared to the actual spectrum (solid line). The x-axis corresponds to MFP for the 
MFP spectrum or the grating period for thermal diffusivity data. The reconstructed spectrum 
agrees well with the actual function provided high spatial frequency data points from (a) are 
discarded. 


grating period, ~ 27t/1 /itttT 1 , are only minimally constrained by the available data, and 
may thus be unreliable. Discarding these points, we perform the reconstruction with the 
subset of data corresponding to spatial frequencies less than ~ 2ir The result is 

given in Fig. 2b, demonstrating quite reasonable agreement with the actual MFP spectrum. 
Thus, the MFP spectrum can be effectively extracted from 2D TG measurements even if 
multidimensional heat conduction occurs in all of the measured data. The reconstruction 
can be made more robust if data with grating period smaller than the optical decay length 
can be obtained so that the heat conduction is primarily one-dimensional for some of the 
measurements. 
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IV. ANISOTROPY OF APPARENT THERMAL CONDUCTIVITY 


A recent experimental work reported the anisotropic failure of Fourier’s law in a time- 
domain thermoreflectance (TDTR) experiment in which heat conduction occurs both in- and 
cross-plane. 28 The authors observed that the in-plane thermal conductivity of Si appeared to 
be nearly a factor of two smaller than the cross-plane value when the pump beam diameter 
was on the order of a few microns, comparable to the phonon MFPs in Si. 

Although the samples in the TDTR experiment consist of a film on top of a substrate 
rather than only a substrate as in TG, we can qualitatively examine this report using our 
multidimensional solution of the BTE to describe heat conduction in the substrate. First, 
consider the form of apparent thermal conductivity expressed in Eq. 5. This equation 
shows that for an isotropic crystal, the apparent thermal conductivity depends only on the 
magnitude of spatial frequency rather than the components along each direction, implying 
that the apparent thermal conductivity must be isotropic in real space. Since Si is a cubic 
crystal that is thermally isotropic, our BTE solution indicates that the apparent thermal 
conductivity must be isotropic even in the quasiballistic regime regardless of the anisotropy 
of the thermal gradients. 

The obvious question is then why an anisotropic thermal conductivity appears to explain 
the data even though the apparent thermal conductivity must be isotropic. A possible 
answer can be found in the similarity of the Fourier’s law solution used to fit the data 
and the exact solution given by the BTE. Recall from Eq. 9 that the temperature decay of 
a 2D grating consists of a multiexponential decay with time constants determined by the 
apparent thermal conductivities that depend only on the magnitude of the spatial frequency. 
On the other hand, the equation giving the Fourier’s law solution for an anisotropic thermal 
conductivity is given by: 

ts_ A' r eX Pb(«^ + <*y(%)t] 

u I (fi+Qm+ey) 4 * 

where a x and a y are the thermal diffusivities along the x and y directions, respectively, 
and are used as fitting parameters to fit the measured thermal decay. Consider the similari¬ 
ties between Eqs. 9 and 11. In Eq. 9, the integral contains a spectrum of thermal diffusivities 
that continuously vary with spatial frequency £. In Eq. 11, there are only two thermal diffu¬ 
sivities but these separately multiply £ x and £ y . Due to the averaging caused by the integral 
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over £ y , it is possible to reproduce a transient thermal decay by appropriately adjusting a x 
and Oi y even though this anisotropic solution is inconsistent with the actual BTE solution. 

As an example, we calculate the thermal decay derived from the BTE using Eq. 9, then 
attempt to fit this decay curve using Eq. 11 with a x and a y as fitting parameters. For 
A = 20 fim and f3 = 1/500 nm _1 , we are able to fit the thermal decay nearly exactly using 
a x = 0.85 abuik and a y /a x = 0.86 as shown in Fig. 3. However, this successful fitting does 
not imply that the apparent thermal conductivity is anisotropic, but rather highlights that 
multiple solutions can explain the same data for this ill-posed problem. This observation 
demonstrates the utility of having the exact analytical solution to the multidimensional BTE 
as presented here. 

Our result shows that anisotropic thermal gradients in bulk isotropic crystals are not 
sufficient to yield apparent thermal conductivities that are anisotropic. It is possible that an 
anisotropic heating profile, in combination with the metal-semiconductor interface present 
in TDTR, could lead to anisotropic apparent thermal conductivity; further study with a 
BTE solution that incorporates the interface is necessary to completely resolve the issue. 


V. SUMMARY 

We have analyzed multidimensional quasiballistic thermal transport in TG using an exact 
analytic solution of the BTE. We find that the thermal decay in multidimensional TG 
consists of a multiexponential decay with the time constants determined by the apparent 
thermal conductivities as a function of spatial frequency. By solving a nonlinear inverse 
problem, we are able to recover these apparent thermal conductivities and thereby obtain 
the MFP spectrum of the thermal phonons without any fitting parameters or knowledge 
of the MFPs. This work extends the MFP spectroscopy method with TG to opaque bulk 
materials. 
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FIG. 3: Calculated temperature versus time from the BTE solution (solid line) and from an 
anisotropic Fourier model (dashed line). The anisotropic Fourier model is able to explain the 
observed thermal decay even though the apparent thermal conductivity is isotropic. 
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